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ABSTRACT 

The density variance - Mach number relation of the turbulent interstellar medium is 
relevant for theoretical models of the star formation rate, efficiency, and the initial 
mass function of stars. Here we use high-resolution hydrodynamical simulations with 
grid resolutions of up to 1024^ cells to model compressible turbulence in a regime sim¬ 
ilar to the observed interstellar medium. We use Fyris Alpha, a shock-capturing code 
employing a high-order Godunov scheme to track large density variations induced by 
shocks. We investigate the robustness of the standard relation between the logarithmic 
density variance (cTg) and the sonic Mach number (M.) of isothermal interstellar tur¬ 
bulence, in the non-isothermal regime. Specifically, we test ideal gases with diatomic 
molecular (7 = 7/5) and monatomic (7 = 5/3) adiabatic indices. A periodic cube of 
gas is stirred with purely solenoidal forcing at low wavenumbers, leading to a fully- 
developed turbulent medium. We find that as the gas heats in adiabatic compressions, 
it evolves along the relationship in the density variance - Mach number plane, but 
deviates significantly from the standard expression for isothermal gases. Our main 
result is a new density variance - Mach number relation that takes the adiabatic in¬ 
dex into account: = In (l -|-provides good fits for bAi < 1. A 

theoretical model based on the Rankine-Hugoniot shock jump conditions is derived, 
tTg = ln{l -I- (7 -I- — l)b'^M^ + 2]}, and provides good hts also for bM > 1. 

We conclude that this new relation for adiabatic turbulence may introduce important 
corrections to the standard relation, if the gas is not isothermal ( 77 ^ 1 ). 

Key words: equation of state - galaxies: ISM - hydrodynamics - ISM: clouds - ISM: 
structure - turbulence 


1 INTRODUCTION 


The interstellar medium (ISM) is a complex, turbulent, 
multi-phase gaseous medium, which permeates the space be¬ 
tween stars in the galactic plane ( Ferriere|2001 1. It is an es¬ 
sential part of the evolutionary cycle in stars, recycling the 
products of nucleosynthesis from dying stars and creating 
the stellar nurseries for a new generation of star formation 
||Mac Low &: Klessen|200^ [Elmegreen fc Scalo|2004[ |McKee| 


fc Ostriker||2007 Krumhoiz||2014 Padoan et aI.||20T4 l. The 

ISM interacts with supernova explosions, protostellar jets, 
winds and outflows, which shape its structure and drive the 
turbulence we observe via atomic and molecular line obser¬ 
vations of the ISM. 


In many simulations that include an ISM (e.g.. Bland- 


[Hawthorn et al.|2007[ [Wagner fc Bicknell|201H[Cooper et al. 
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2009 Fischera & Dopita 20051, a statistical construction 


with isotropic properties related to turbulent statistics have 
been used, as a proxy for the turbulent ISM. Causal models 
of the turbulent ISM will drastically increase the accuracy 
of these models, but this first involves an in-depth study of 
the statistics and evolution of fully-developed turbulence. 

In purely isothermal gas, the probability density func¬ 
tion (PDF) of the gas densities may be approximated by a 
lognormal distribution ( |Vazquez-Semadeni||1994| ), which in 
log-space has the form 


Pln(s) = 




: exp 


1 - s? 

2 0-2 


( 1 ) 


where s = ln(p/p), s and are the mean and variance of 
the logarithm of the density p, scaled by the mean density, 
p. The logarithmic density variance is a function of the root- 
mean-squared (rms) sonic Mach number (AI), and is given 
by 


= In (1 -I- ■ 


( 2 ) 
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The coefficient b is known as the turbulence driving param¬ 
eter and depends on the mode mixture induced by the tur¬ 
bulent forcing mechanism (Federrath et al. 2008|). Purely 


solenoidal (divergence-free) driving leads to b — 1/3, while 
purely compressive (curl-free) driving corresponds to 6 = 1 . 
Equation jiL has been studied extensively for isother- 


mal gases (Padoan et al.|; 

997| Passot & Vazquez-Semadeni 

1998 

IKritsuk et al.|20071 

Beetz et al.|2008 I’ederrath et al. 

2008 

Price et al. 201 1| | 

Burkhart & Lazarian| 2012 Seon 

2012 

|Konstandin et al.|2012l, with investigation into differ- 

ent simulation techniques (Price & Federrath|2010 1 and stir- 

ring methods (Federrath et al.|2008 2010). It has also been 


studied by employing a heating and cooling curve (IWada & 

Norman|2001 

Kritsuk & Norman|2002| 

Audit & Hennebelle 

2001" mToT" 

ennebelle & Audit|f20071 

Seifried et al. 20111 

Gazol & Kim 

20131. Recently an investigation has been done 


in the magnetohydrodynamic (MHD) regime (Molina et al. 
2012 1 , and on polytropic gases (Federrath & Banerjee 20151. 


In our work, we investigate the robustness of this well- 
established density variance - Mach number relation. Equa¬ 
tion in the non-isothermal regime, specifically in ideal 
gases with diatomic molecular (7 = 7/5) and monatomic 
(7 = 5/3) adiabatic indices. 

This is relevant because theoretical models of the 


star formation rate (Krumholz & McKee 


2005 


Padoan & 


Nordlun^|201H |Hennebelle fc Chabrier||2011[ |Federrath~^ 


Klessen 


_ 2012 1, the star formation law (Federrath 2013bl, 

the star formation efficiency ( Elmegreen| 2008 1 , and the ini¬ 
tial mass function of stars ([Hennebelle fc Chabrier||2008[ 


Hopkins||2013a Chabrier et al.||2014 l heavily rely on Equa- 

tion ([ 2 |). 

Section [^summaries our simulation and analysis meth¬ 
ods, Section^ first presents results for the isothermal case, 
in order to make contact with previous studies and to verify 
our analysis techniques. Then we present a numerical reso¬ 
lution study to determine the minimum resolution required 
in order to measure the density variance - Mach number 
relation in simulations with 7 > 1 and present our main re¬ 
sults for adiabatic indices 7 = 7/5 and 5/3. We provide a 
theoretical model for the crJ(Af) relation in Section]^ and 
discuss the discrepancies that we find compared to the stan¬ 
dard Equation (§. Section summarises our conclusions. 


2 SIMULATION AND ANALYSIS METHODS 


To simulate the turbulent ISM we use the high-resolution. 


shock-capturing code Fyris Alpha (Sutherland 2010 1 to solve 


the equations of compressible hydrodynamics across a three- 
dimensional, periodic domain with side length L = 1 , initial 
uniform density p = 1 , pressure of 1/2 (Cg = 7 / 2 ), and 
zero initial velocities. Unlike previous studies, the goal here 
is to test the density and velocity statistics of purely adia¬ 
batic turbulence with an ideal gas EOS, rather than a purely 
isothermal or polytropic EOS, and simpler than employing a 
cooling curve or running chemo-hydrodynamical simulations 


(Glover et al. 2010). Ultimately, simulations with multiple 


species including all relevant chemical reactions, as well as 
radiative heating and cooling would be the most realistic, 
but their complexity might not allow us to reduce the results 
to some simple rules of thumb that can be used in practical 
applications. We thus simplify the problem significantly by 


studying the turbulence in purely adiabatic, ideal gases with 
the aim of extracting results that might be applicable to a 
wider range of cases, including terrestrial experiments and 
atmospheric turbulence, in addition to the ISM. Tablej^lists 
the key parameters of all our adiabatic turbulence simula¬ 
tions. 


2.1 Ideal gas equation of state 


The ideal gas equation of state relates the pressure P, den¬ 
sity p and temperature T, and is given by 

PV P N 

—- = Nku or — = —kuT or P = nkBT, (3) 
T pm 

with the volume U, the Boltzmann constant fca, the particle 
mass m, the total number of particles A, and the particle 
number density n = N/V. The ratio of the specific heat 
capacities at constant pressure and constant volume defines 
the adiabatic index. 


_ Nt 

cv 



(4) 


where / denotes the number of degrees of freedom. For 
monatomic gas, / = 3 and 7 = 5/3, while for diatomic 
molecular gas, / = 5 and 7 = 7/5, because diatomic 
molecules have two rotational degrees of freedom in addition 
to the three translational degrees of freedom. Note that at 
typical molecular cloud temperatures (about 10-100 K), os¬ 
cillatory degrees of freedom cannot be excited by collisions, 
which is why—although theoretically present—they do not 
contribute to increase / in such cases. The specific internal 
energy of an ideal gas, u = is only determined by 

its temperature. Inserting this equation into Equations § 
and (|^, leads to 

P{p,T) = {j - l)pu{T) (5) 


expressed via the adiabatic index 7 , which serves as the 
equation of state. In order to determine the statistics of tur¬ 
bulence in this adiabatic regime, we use isothermal, diatomic 
molecular and monatomic equations of state ( 7 —>-1 ,7 = 7/5 
and 5/3 respectively). 

Note that in order to model isothermal gases, 7 is often 
set close to unity (e.g., 7 = 1 . 0001 ), as if the gas had an 
extremely large number of degrees of freedom / —^ 00 . This 
trick produces a gas that approximately stays at constant 
temperature, because any excess heat from dissipation (e.g., 
by shocks) is absorbed in such a big internal energy reservoir 
that the temperature of the gas does not notably change. 


2.2 Driving of turbulence, time evolution, and 
definition of the Mach number 

The driving of turbulence in the gas is performed by stir¬ 
ring with random purely solenoidal (divergence-free) forcing 
at low wavenumbers for the duration of the simulation. All 
wavenumbers k in the range 1 ^ fc/(27r/Z/) ^ 3 were driven. 
The driving pattern is evolved with an Ornstein-Uhlenbeck 
process similar to the methods exp lained in|Eswaran fc Pope| 
(1988), Schmidt et al. (2009), and Federrath et al. ( 2010| ). 

From the work done by Federrath et al.| ( 2008 2010[ ), we 
expect a proportionality constant of b ~ 1/3 in Equation it 
for solenoidally driven isothermal gas and we will test that in 
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Table 1. Simulation parameters. Notes. Column 1: simulation name. Column 2: grid resolution. Column 3: adiabatic exponent 7 in 
Equation Column 4: dimensionless driving amplitude of the turbulence. Column 5: Resulting time-averaged velocity dispersion in 
code units in the regime of fully developed turbulence. Column 6 : turbulent box crossing time: tcross = L/cr^ in code units. 


Simulation name 

7V3 

'^res 

7 

A 

G y 

Across 

(01) AD-TURB-256-A200-G1 

256^ 

1.0001 

200 

2.19 ±0.14 

0.46 ±0.03 

(02) AD-TURB-256-A400-G1 

2563 

1.0001 

400 

3.38 ±0.13 

0.30 ±0.01 

(03) AD-TURB-256-A800-G1 

2563 

1.0001 

800 

5.45 ±0.29 

0.18 ±0.01 

(04) AD-TURB-256-A1600-G1 

2563 

1.0001 

1600 

8.62 ±0.61 

0.12 ± 0.01 

(05) AD-TURB-256-A100-G7/5 

2563 

7/5 

100 

1.53 ±0.04 

0.65 ± 0.02 

(06) AD-TURB-256-A200-G7/5 

2563 

7/5 

200 

2.47 ± 0.12 

0.40 ±0.02 

(07) AD-TURB-256-A400-G7/5 

2563 

7/5 

400 

4.00 ±0.17 

0.25 ± 0.01 

(08) AD-TURB-1024-A200-G5/3 

10243 

5/3 

200 

2.51 ± 0.14 

0.40 ± 0.02 

(09) AD-TURB-512-A200-G5/3 

5123 

5/3 

200 

2.50 ±0.13 

0.40 ± 0.02 

(10) AD-TURB-256-A100-G5/3 

2563 

5/3 

100 

1.55 ±0.05 

0.65 ± 0.02 

(11) AD-TURB-256-A200-G5/3 

2563 

5/3 

200 

2.55 ±0.14 

0.39 ±0.02 

(12) AD-TURB-256-A400-G5/3 

2563 

5/3 

400 

4.07 ± 0.24 

0.25 ± 0.01 

(13) AD-TURB-128-A200-G5/3 

1283 

5/3 

200 

2.47 ± 0.13 

0.40 ± 0.02 

(14) AD-TURB-64-A200-G5/3 

643 

5/3 

200 

2.45 ±0.14 

0.41 ± 0.02 



Figure 1. The velocity dispersion (t„, as a function of simulation 
time for driving amplitudes A = 100, 200 and 400, and for fixed 
adiabatic 7 = 7/5. The times for the onset of turbulence in each 
case are shown as vertical dashed lines and are approximated 
with the box crossing time tcross = T/(t„, where L is the linear 
size of the computational domain. The box crossing time for each 
simulation is listed in Table [T] 


both isothermal and adiabatic gases. The rms Mach number 
of the gas is modified by varying the stirring amplitude A 
of the driving force, allowing each 7 to be tested at a range 
of Mach numbers. 

All simulations are run for several turbulent box cross¬ 
ing times to test the density variance - Mach number re¬ 
lation in the regimes of transient as well as fully-developed 
turbulence. The turbulent crossing time is defined as tcross = 
LjcTv, where (t„ is the asymptotic velocity dispersion. The 
velocity dispersion and hence the crossing time depend on 
the driving amplitude A. We show an example of this de¬ 
pendence in Fig. [2 We assume that the turbulence becomes 
fully developed after one crossing time in each simulation, 


indicated by the dashed vertical lines in Fig.[^ At this point 
the gas properties no longer vary drastically but change 
smoothly, indicating a statistically stable configuration. This 
allows us to distinguish regimes of transient {t < tcross) and 
fully-developed (t > tcross) turbulence. 

Given the velocity dispersion and sound speed Ca = 
{dP/dpy^^, we have different Mach numbers A4 = avjcs, 
depending on the driving amplitude A and depending on 
the value of adiabatic 7 . This is because the sound speed 
depends on the derivative of the EOS, Equation fit- It fur¬ 
thermore depends on the simulation time, because the in¬ 
ternal energy and thus the mean temperature of the gas 
keeps increasing during the course of the adiabatic simula¬ 
tions with 7 = 7/5 and 7 = 5/3. This is in stark contrast 
to the isothermal and polytropic simulations performed in 


previous studies ( 

Padoan et al. 

1997 

Passot & Vazquez- 

Semadeni|19981|Federrath et al.|2008 Price et al.|2011||Kon- 

standin et al.||2012 Federrath & Banerjee||2015 

1 , where the 


sound speed did not change systematically, after the turbu¬ 
lence was fully developed. Here, however, the total energy 
is conserved, which means that all dissipated energy is con¬ 
servatively added to the internal energy. Thus, the injected 
energy from the driving is converted into internal energy 
and heats the gas continuously, leading to an ever increasing 
average sound speed and to a continuously decreasing rms 
Mach number. We thus use instantaneous measurements of 
M and a^, in the following to determine the density variance 
- Mach number relation in adiabatic gases. 


2.3 Measuring the density variance 

The density variance of the gas is calculated using method 
4 in Section 2.3 of Price et al. (2011 1 , but instead of fitting 
a lognormal distribution, we fit the more appropriate |Hop^ 
|kins| | |2013b[ ) distribution. The advantage of the Hopkins ht 
is that it takes turbulent intermittency effects into account 
and provides excellent hts to the density PDFs over a wide 
range of physical parameters, including different Mach num¬ 
bers, driving amplitudes and mixtures |Federrath|2013al, as 
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Mach number 


Figure 2. The density variance — Mach number relation for 
isothermal turbulence (approximated by setting 7 = 1.0001). In 
order to reach a wide range of Mach numbers to test the re¬ 
lation, we use stirring amplitudes between A = 200 and 1600 
(models 1-4 in Table [l)| , leading to Mach numbers between 3 and 
12. We fit Equation | [^ to the four points and find the proportion¬ 
ality constant b = 0.37di 0.10 with a goodness of fit of = 0.99. 
The fit is shown as a solid line, while cases with b = 0.3 and 
b = 0.4 are shown as dashed lines for comparison. Our best fit is 
consistent with the expectation value 6 ~ 1/3 for purely solenoidal 
driving ( [Federrath et al.|2008|[^10| |. 


well as magnetic field strengths and variations in the poly¬ 
tropic exponent for simulations that employ a polytropic 
EOS (Federrath & Banerjee 20151. The Hopkins (2013bI 
density PDF is defined as 


Phk(s) = h (2^/Xuj{s)^ exp [- (A + a;(s))] ^ 


X 


9^ uj{s) 


X = a7(20"), Ljis) = A/(l + 9)- s/9 {u ^ 0), ( 6 ) 


where li{x) is the modified Bessel function of the first kind. 
Equation (|^ is motivated and explained in detail in |Hop- 


kins ( 2013b||. It contains two parameters: 1) the volume- 


weighted standard deviation of logarithmic density fluctua¬ 
tions (Ts, and 2) the intermittency parameter 9. In the zero- 
intermittency limit (0 —>■ 0), Equation § becomes the log¬ 
normal distribution from Equation 0 , Phk —>■ Pln- 

In order to measure the density variance we fit our 
simulation density PDFs in a restricted range around the 
mean (from s —Sus.mom to s-b Sus.mom, where (T 3 ,mom is the 
second moment of the density distribution, directly com¬ 
puted by summation over all simulation data points) with 
Equation 0 and determine the best-fit parameter as- In 
agreement with the conclusions drawn in [Price et al. (2011 1 
and Hopkins (2013bI, we find that the fitted is the same 
within a few percent as (computed by summation 

over all simulation grid cells). 


3 RESULTS 

3.1 Isothermal comparison 

The isothermal gas relation has been studied extensively and 
is therefore a good comparison test for our hydrodynamics 
code, setup and post-processing methods to determine Us 
(see 0 - For purely solenoidal forcing of the turbulence, we 
expect a proportionality value of 6 ~ 1/3 in Equation 0 
( [Federrath et al^|2008| |2010p . Four simulations were per¬ 
formed at a grid resolution of 256® with 7 = 1.0001 to pre¬ 
vent non-isothermal effects brought on by high stirring am¬ 
plitudes. The four time-averaged points lie between b — 0.3 
and 0.4, with a fit of 6 = 0.37 ± 0.10 and goodness-of-fit pa¬ 
rameter R® = 0.99 (Fig. § . Our measurement of b spans the 
expected value, thus our methods produce reasonable results 
for isothermal turbulence. Now that we have established that 
our simulation and analysis techniques reproduce previous 
results for isothermal turbulence, we can now move on to 
study non-isothermal turbulence in the adiabatic regime. 


3.2 Resolution study 


We now test the resolution requirements of density variance 
and Mach number measurements by performing a series of 
identical simulations at increasing resolutions of 64®, 128®, 
256®, 512® and 1024® grid cells for 7 = 5/3. Fig.|^shows the 
density variance - Mach number relation for each simula¬ 
tion and simulation time. We see that first, the density vari¬ 
ance and Mach number increase and reach a maximum after 
about tcross (the lower part of the correlation). In this first 
part of the evolution, the kinetic energy of the gas increases 
due to the driving until the kinetic energy power spectrum 
is established ( Schmidt et al.|2009 l. After the kinetic energy 
and rms velocity have reached a saturated state (see Fig.[^, 
only the sound speed keeps increasing monotonically, be¬ 
cause the dissipated energy heats the gas. This leads to a 
continuously decreasing Mach number and density variance 
in the regime of fully developed turbulence (the upper part 
of the correlation). 

The higher the numerical grid resolution, the greater 
the maximum of uj for Wes < 256, which is seen to con¬ 
verge for Wes ^ 256. A zoom within the focus area of the 
density variance - Mach number relation is shown in the 
right-hand panel of Fig. [^ These points are independent of 
the initial jump, but are still seen to converge at increased 
resolution. For resolutions equal to and above 256®, the val¬ 
ues are almost identical. Therefore the statistics of density 
variance and rms Mach number may be approximated well 
by numerical resolutions ^ 256® cells. This is consistent with 
the resolution requirements est ablished in [Kitsionas et al.| 
(20091, Federrath et al. ( [2010[ ), [Kritsuk et al. (20111, and 
Federrath (2013a I. 


Looking closely at Fig. [^ we see that the gas evolves 
along a curve in the density variance - Mach number plane. 
This is due to the continuous heating of the gas via the 
turbulent driving, which lowers the Mach number continu¬ 
ously. The evolution of this curve seems to correlate some¬ 
what with Equation 0 : but is steeper than the theoretical 
prediction for isothermal turbulence. The behaviour of this 
curve is quantified in detail in the following subsections, §3.3| 


and (3.4 for 7 = 7/5 and 7 = 5/3, respectively. 
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Figure 3. Left panel: Evolution of simulations with increasing resolution in density variance - rms Mach number space. The dashed 
lines represent functions of Equation with b = 0.3 and 0.4, for comparison. Initially the rms Mach number increases substantially, 
then decreases smoothly once the turbulence becomes fully developed (at about tcross), as a result of the continuously increasing internal 
energy, temperature and sound speed for our ideal gas EOS, Equation Right panel: A magnified region of the left panel, displaying 
the convergence of statistics. Simulations with ^ 256^ grid cells are representative of the converged system. 



Figure 4. Evolution of simulations with 7 = 7/5 and A = 100 
(red, small arrows), 200 (purple, normal arrows) and 400 (blue, 
large arrows). The arrows indicate the direction of time evolution. 

3.3 Diatomic molecular gas: 7 = 7/5 

In the first case we look at a diatomic equation of state, 
i.e., 7 = 7/5. We perform three separate simulations with 
different stirring amplitudes A, and plot their evolutionary 
curves in Fig. The best fit to Equation it from the com¬ 
bined points of all three simulations is b — 0.37 ± 0.02 with 
a goodness of fit parameter of = 0.90. Note that the in¬ 
crease in sound speed due to gas heating is slow compared 
to the time it takes to establish a statistically steady state, 
which is shown in Fig. where we see that the velocity 
dispersion is fully developed after one crossing time. How¬ 
ever, the continuous heating in the fully developed regime of 
turbulence leads to a continuously increasing sound speed, 
the consequences of which are discussed in more detail in 
Section [4] 


As in Fig.[^ we see in Fig.fflthat simulations with 7 > 1 
produce a somewhat steeper (jj(A4) relation compared to 
the isothermal relation (cf. Fig. as time progresses and 
both (tJ and M decrease due to the continuous heating of 
the gas. The different times of each simulation are connected 
by a line with arrows in Fig. indicating increasing time. 
Therefore a new functional fit might be more appropriate 
to describe the behaviour in our non-isothermal, adiabatic 
turbulence simulations. 

The simplest modihcation to the existing model func¬ 
tion, Equation is to allow for variations in the exponent 
on the Mach number. Our data in Fig. indicate that the 
exponent is somewhat higher than the standard depen¬ 
dence from Equation §. Thus we use the following new fit 
function to determine the exponent a: 

=ln{l + b'^M"). (7) 

We do not necessarily expect that the coefficient b' in this 
new relation is the same as b in Equation ([^, but we will 
test that below. The new function is fitted to the data from 
each of the three simulations with different driving ampli¬ 
tude and to the combined set of data points. We determine 
the goodness of £t parameter for the fits to Equations ([^ 
and Q and compare them. The results are summarised in 
Tabl^ 

Table[2]shows that for the individual simulations as well 
as for the combined data set, the modified power-law func¬ 
tion from Equation 0 provides the better fit to the data as 
quantified by the goodness of fit R^. The coefficient value 
b' — 0.31 ± 0.04 is smaller than b = 0.37± 0.02, but they are 
formally consistent with representing the same value, and 
consistent with the fe-value obtained in our isothermal cal¬ 
culations in §3.1| In fact, our new fit gives a value that is 
in agreement with the theoretical expectation for the turbu¬ 
lent driving, namely 6 ~ 1/3 for purely solenoidal driving as 
applied here. 

The exponent a, which is Hxed to q = 2 in Equation (2 1 , 
but allowed to vary in our new fit function, Equation (71, 
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Table 2. Statistical fit parameters for different functions crJ(A4), for 7 = 7/5. Notes. R? denotes the goodness-of-fit parameter. A value 
of = 1 indicates a perfect fit to the given data. 


Fit function 

(T^ = ln(l - 1 - b'^AF) 

0-2 = ln{l - 1 - b' 2 At“) 


Parameter 

b 


b' 

a 


A = 100 

0.35 ±0.02 

0.90 

0.32 ±0.05 

2.7 ± 0.8 

0.96 

A = 200 

0.38 ±0.04 

0.90 

0.31 ± 0.07 

2.9 ±0.8 

0.99 

A = 400 

0.37 ± 0.04 

0.92 

0.30 ±0.07 

2.8 ±0.8 

1.0 

All data 

0.37 ± 0.02 

0.90 

0.31 ± 0.04 

2.8 ±0.4 

0.99 



Figure 5. Evolution of simulations with 7 = 5/3 and A = 100 
(red, small arrows), 200 (purple, normal arrows) and 400 (blue, 
large arrows). The arrows indicate the direction of time evolution. 

clearly shows that an almost cubic dependence on M pro¬ 
vides a better fit to the data. We find a best-fit value of 
a = 2.8 ± 0.4 in our simulations with 7 = 7/5, leading to 
a new form of the density variance - Mach number relation 
for 7 = 7/5 gas, 

as = In [1 + (0.31 ± 0.04)^ • x(2'8±0'4)j ^ (g) 

This result presents an interesting question: is the den¬ 
sity variance - rms Mach number relation of non-isothermal 
adiabatic turbulence no longer a quadratic relation, com¬ 
pared to the well-studied isothermal case? We will now ex¬ 
plore whether the same/similar holds for 7 = 5/3 and then 
address this questions in the discussion of Section]^ by com¬ 
paring to a theoretical model of the ( 7 g(A 4 ) relation. 

3.4 Monatomic gas: 7 = 5/3 

In the second case we look at a monatomic EOS, Equa¬ 
tion with 7 = 5/3. As in the diatomic case we performed 
three separate simulations with different driving amplitude 
A, and plot their evolutionary curves in Fig. (again with 
arrows indicating the continuous time evolution to smaller 
and smaller a^ and A4). A fit with the standard relation. 
Equation § , to the combined data set gives 6 = 0.36 ±0.02 
with a goodness of fit = 0.90. 

Given the same effect occurs as for 7 = 7/5, we fit our 
new power-law function. Equation Q, to data from each of 



Figure 6. Exponent o in the density variance - Mach number 
relation for different values of the adiabatic index 7 . The data 
points show our simulation measurements and the dotted line is 
a linear fit with a = (57 -I- l)/3. 

the three simulations and to the combined set of data points, 
and compare the goodness of fit values with those from the 
fit to Equation ([^. The results are summarised in Table 
Once again we see a significantly better fit to the power 
law with exponent a > 2. We find that the driving coefficient 
b' ~ b, as for the 7 = 7/5 case, indicating that the physics of 
the driving is indeed contained in the value of the b param¬ 
eter, while the fact that we deal with non-isothermal tur¬ 
bulence is reflected in a steeper power-law exponent a > 2 , 
compared to the isothermal case. All three simulations fit 
exponents very close to cubic {a ~ 3), with the combined 
data for 7 = 5/3 fitting a power law of the form 

crj = In [1 ± (0.32 ± 0.03)^ • x (3 0 ±o. 5 )j ^ 


3.5 Summary of the results 

In summary, both the 7 = 7/5 and 7 = 5/3 cases yield 
turbulent driving coefficients 6 ~ 1/3 that are all consistent 
with the theoretical expectation for purely solenoidal driving 
of the turbulence, and consistent with the 6 -values found for 
isothermal turbulence (7 —> 1 ). 

The exponent a of the aRA4) relation, however, is 
significantly steeper with a ~ 2.8 ± 0.4 for 7 = 7/5 and 
a = 3.0 ± 0.5 for 7 = 5/3, compared to the isothermal case, 
where a = 2 provides the best fit to the data. Thus, we see 
that the exponent a in the density variance - Mach number 
relation is 7 -dependent. 
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Table 3. Statistical fit parameters for different functions ag(A4), for 7 = 5/3. Notes. R? denotes the goodness-of-fit parameter. A value 
of = 1 indicates a perfect fit to the given data. 


Fit function 

0-2 = ln(l -1- b’^JVp) 

0-2 = in(i 


Parameter 

b 

A2 

b' 

a 

R2 

A = 100 

0.34 ± 0.03 

0.92 

0.32 ±0.05 

2.8 ± 0.9 

0.99 

A = 200 

0.38 ±0.04 

0.91 

0.33 ±0.06 

2.9 ± 0.8 

1.0 

A = 400 

0.36 ±0.04 

0.89 

0.32 ±0.06 

3.0 ± 0.9 

1.0 

All data 

0.36 ±0.02 

0.90 

0.32 ±0.03 

3.0 ± 0.5 

0.99 


In order to provide a heuristic relation that describes 
the behaviour of in adiabatic gases, we fit a linear 

function to the dependence of a on 7 in Fig. It is rea¬ 
sonable that a will continuously increase with 7 . Thus, we 
chose the simplest function to approximate our data from 
7 = 1 to 7 = 5/3, i.e., a linear interpolation. The result is 
a = ( 57 -I- l)/3. The actual dependence might be somewhat 
different, but to determine a better fit would require us to 
measure a for a range of 7 values, in much smaller steps A 7 . 
This is beyond the scope of the paper, but we can already 
provide a new improved functional form of the density vari¬ 
ance - Mach number relation that takes the adiabatic index 
7 into account. The best-fit function is given by 

+ ( 10 ) 

which is the central result of the paper. Equation ( |10| | nat¬ 
urally simplifies to the well-studied isothermal case (7 —>■ 1 ) 
given by Equation but also approximately covers cases 
with 7 > 1, up to 7 = 5/3. 


4 DISCUSSION 


In §3.3| - g375l we found that the density variance - Mach 
number relation for adiabatic gases with 7 = 7/5 and 5/3 
respectively, deviates significantly from the isothermal case. 
Equation (§. We quantified the discrepancy by fitting an 
alternative function. Equation 0 , to the data, with the 
power-law exponent a as a free fit parameter. We find that 
the power law provides excellent fits, with power-law expo¬ 
nents increasing with 7 from a = 2 for the isothermal case 
(7 —>■ 1) to a ~ 3 for 7 = 5/3. A heuristic function was 
obtained to provide a new crg(Af) relation that takes the 
dependence on 7 into account, given by Equation (lOl. 


We now compare this results to a recent theoretical 
model for the density variance - Mach number relation, in 
order to explain the differences of our adiabatic case to the 
isothermal and polytropic cases. The detailed derivation of 


the relation can be found in Molina et al. (20121 and Fed 
[errath fc Banerjee] ( |2015[ ), where this relation has been ex¬ 
plored for magnetized isothermal and polytropic gases re¬ 
spectively, with the latter representing a special case, where 
the pressure and temperature of the gas are both uniquely 
related to the density via 


P{P) 


Tip) 


( 11 ) 


We emphasize that this is different from the adiabatic EOS, 
Equation Q, used here, in that the pressure depends on 
both density and temperature, P{p, T). Thus, for any given 


value of density, the pressure can vary depending on the 
temperature, while a polytropic EOS will give only a single 
value of P for a given input p. 

The basic idea of the theoretical model is to relate 
the density jump in a single shock to the ensemble of 
shocks/compressions in a turbulent medium. For that pur¬ 
pose, [Padoa^^^^ordlu^ (20111 and Molina et al. (20121 
applied the equations of mass, momentum and energy con¬ 
servation across a shock, 


Povo = pv, 

povo -I- Po = pv^ + P, 

I2. I2, ,P 

t:Vo+uo-\ - =-V -hu-l-, 

2 po 2 p 


( 12 ) 

(13) 

(14) 


to derive the density contrast p/po between the pre-shock 
gas (denoted with index 0 and on the left-hand side of the 
equations) and the post-shock gas (no index; right-hand side 
of the shock jump equations). The equation of state. Equa¬ 
tion 0 . enters through the pressure P and specific internal 
energy u in these expressions. Combining these equations 
leads to the well-known Rankine-Hugoniot shock jump con¬ 
ditions as the solution for the density jump across the shock 


( Rankine|1870 Hugonio't||1887| Shull fc Draine[|1987 l 

P _ Vo _ ( 7 - 1 - 

Po 


Vo 

V 


( 7 - 1)62X2 -t 2 ’ 


(15) 


Note that we have already introduced the geometrical 6 - 
parameter, because these shock jump conditions only apply 
to the plane-parallel component of the shock, parametrized 
by the parallel component of the sonic Mach number 


vo/csfi = bA4 (Molina et al. 2012 Federrath & Banerjee 
20151. Following the detailed derivation in Molina et al. 


(20121, Equation (15 \ just needs to be inserted into the gen¬ 


eral expression for the ensemble of such shocks. 


crj = In ( 1 -I- -((- 

po 


(16) 


which leads to the density variance - Mach number relation, 

i'y + l)b^M^ 


CTg = In ( 1 -I- 


(7 - 1)62X2-h 2 


(17) 


We immediately see that this new 7 -dependent density vari¬ 
ance - Mach number relation reduces to the isothermal case. 
Equation 0 , if we set 7 = 1. In the adiabatic case, how¬ 
ever, 7 > 1 , which leads to our theoretical prediction as a 


function of 7 , given by Equation (17l. 


Fig. [7] shows the theoretical prediction given by Equa¬ 
tion flTf for different values of 7 together with the isother¬ 
mal solution and together with our simulation data for 
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Figure 7. Combined density variance - Mach number relation plot, showing all our simulation data with 7 = 1.0001 (red crosses), 
7 = 7/5 (green boxes), and 7 = 5/3 (blue diamonds). The dashed lines are our theoretical prediction given by Equation with the 
respective values of 7 (in the same colour as the simulation data and labelled on each curve). 


7 = 7/5 and 7 = 5/3. We see that the new relation quali¬ 
tatively follows the trend of a slightly steeper rise with in¬ 
creasing 7 for low Mach number. It also predicts that at 
high Mach number, the density variance saturates at lower 
cTj for increasing 7 . This is reasonable, because the density 
jump across shocks reduces significantly with increasing 7 
as derived in Equation (151. This regime needs to be tested 
in follow-up simulations that reach higher Mach numbers. 
However, the problem is that the adiabatic heating increases 
with increasing Mach, snch that it quickly counteracts the 
effect of an increased driving amplitude. 


Despite these reasonable qualitative trends produced by 
our new theoretical relation. Equation ( |17| , we see that the 
actual simulation data with 7 > 1 still follow a somewhat 
steeper curve at low Mach number, bM < 1 in the 
plane. We speculate that this discrepancy arises, because the 
theoretical model does not contain any information about 
the temporal evolution of the gas, in particular about its 
temperature changes along the evolutionary curve. 


However, we can qualitatively argue that any shock will 
immediately experience the temperature and pressure in¬ 
crease associated with the adiabatic compression. This will 
reduce the density jump significantly, such that the density 
variance will be smaller with increasing 7 almost immedi¬ 
ately when these shocks are about to form, (e.g., the theo¬ 
retical limit for 7 = 5/3 is p/po = 4). Thus, shocks in high -7 
gas will be significantly reduced and so will be the statisti¬ 
cal variance of density fluctuations, uj. The important point 


here is that this process almost instantaneously reduces the 
density variance. 

At the same time, each shock dissipates energy, locally 
increasing the temperature and internal energy of the ideal 
adiabatic gas. This increases the sound speed, but at first 
only locally in the shocks, which leads to a decreasing pre¬ 
shock sonic Mach number over time, resulting in the time 
evolution (shown as arrows) in Figs. and We can thus 
qualitatively understand the time dependence and result¬ 
ing CTg(A4) relations for 7 > 1. The aRAi) relation is 
steeper, because uj responds almost instantaneously to the 
local pressure and temperature increase in shocks, while the 
Mach number reduction is delayed, because the sound speed 
increases only in the post-shock gas, while our theoretical 
Equation is based on the large-scale, volume-weighted 
pre-shock Mach number. 

In order to substantiate our findings, we show density, 
temperature, pressure, and entropy slices of our highest- 
resolution simulation with 1024^ grid cells and adiabatic 
7 = 5/3 in Fig. We see two important points. First, the 
pressure and temperature are not unique functions of the 
density, but for a given density, the gas can have a range 
of temperatures and pressures, as implied by Equation ©• 
This is substantiated by the entropy slice shown in the bot¬ 
tom right panel of Figure which is not uniform, demon¬ 
strating that the gas is not barotropic, but that the pressure 
depends on density and temperature. There is clearly vis¬ 
cous heating, which is primarily due to shocks in the bAi > 1 
regime, while turbulent dissipation (eddy viscosity) becomes 
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Figure 8. Slices of the normalized density (top left), temperature (top right), pressure (bottom left), and specific entropy (bottom right) 
at f = tcross for 7 = 5/3, A = 200, and a numerical resolution of 1024®. It is evident that gas at a given density can have a wide range 
of temperatures and pressures, unlike a polytropic EOS where T and P are unique functions of the density only. 


a more important heating mechanism when bAi < 1. Quan¬ 
tifying both contributions is beyond the scope of this pa¬ 
per. Second, the adiabatic heating primarily occurs in the 
post-shock gas. The rise of the internal energy does not im¬ 
mediately reduce the global post-shock Mach number, but 
slowly diffuses to large scales, before it affects A4, leading 
to the steeper-than-isothermal as(A4) relations we found for 
7 > 1. 


Finally, Fig. shows density-temperature correlation 
probability density functions (PDFs). It is evident that for 
any given density, there is a wide range of temperatures 
and that the heating indeed primarily occurs in the densest 
gas, i.e., in the post-shock regions. Also note the continuous 
increase in the overall temperature of the gas between t = 
tcrosa (left-hand panel) and t = 2 tcroas (right-hand panel), 
which leads to the slowly decreasing A4 over time. 


5 SUMMARY AND CONCLUSION 


We performed hydrodynamical simulations of supersonic 
and subsonic turbulence, employing an ideal equation of 
state (EOS) with adiabatic indices 7 = 1.0001 (nearly 
isothermal), 7 = 7/5 (diatomic molecular gas), and 7 = 5/3 
(monatomic gas). Sectionprovided a detailed analysis of 
the density variance - Mach number relation, ctJ/AI), which 
is a key ingredient for theoretical models of the star for¬ 
mation rate and the initial mass function. Unlike previous 
studies of purely isothermal and polytropic turbulence, we 
find that an ideal gas EOS leads to a steeper dependence of 
the density variance on the rms sonic Mach number A4. 
We find a new combined approximate relation of the form 


given by Equation (10 1 for low Mach numbers, bA4 < 1, 


which reduces to the well known isothermal solution for the 
special case 7 —>■ 1, but also covers cases 7 > 1. We argue 
that the steeper-than-isothermal dependence for bM < 1 is 
a result of the local heating of the gas in post-shock regions, 
with the global sonic pre-shock Mach number in the relation 
being affected later in the evolution. This is because the tur- 
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Figure 9. Density—temperature correlation PDFs for our adiabatic turbulence simulations with 7 = 5/3, A = 200 and //res = 1024. The 
left-hand panel shows the results at t = tcross, while the right-hand panel is for t = 2tcross. The distributions are spread by more than 
an order of magnitude for typical gas densities, around the average correlations (shown as white lines) with T ~ The continuous 

heating of the gas indicated by the overall rise i n te mperature at t = 2 tcross compared to t = tcross primarily occurs in the post-shock 
gas, while the Mach number entering Equation l|l7| only applies to the pre-shock gas. This produces the steeper dependence of cr^ on 
A4 that we find in our central result. Equation (|10|. 


bulent driving keeps increasing the internal energy reservoir, 
leading to an ever increasing global sound speed in adiabatic 
gases. This is in stark contrast to isothermal and polytropic 
gases, where the sonic Mach number reaches a statistical 
steady state rather than continuously decreasing. 

We derived a theoretical model, Equation (171, for the 
ctRAA) relation in Section]^ which is based on the Rankine- 
Hugoniot shock jump conditions and provides reasonable hts 
to all our data. It furthermore predicts a saturation of crj 
for bA4 >> 1 , which is not yet in reach by numerical sim¬ 
ulations. Such a saturation is reasonable for 7 > 1, given 
the fact that adiabatic shocks always have a finite jump 
in density, while isothermal shocks can theoretically have 
an inhnitely large jump in density across the shock. Both 
Equation (10 1 for bA4 < 1 and Equation (171 for bA4 > 1 
naturally simplify to the standard isothermal relation, Equa¬ 
tion it for 7 = 1 . 

We conclude that changes in the adiabatic exponent 7 
can introduce important modifications in the density vari¬ 
ance - Mach number relation and we provide an approxima¬ 
tion of that behaviour in Equation ( |10[ ). However, we empha¬ 
size that the real ISM is a mixture of atomic and molecular 
phases and that the effective EOS is determined by a com¬ 
plex balance of heating and cooling processes, which in turn 
depend on the chemical evolution and exposure to interstel¬ 
lar and local radiation helds (e.g., from massive stars). Thus, 
our systematic study sheds some light on the dependence of 
turbulent density fluctuations on the thermodynamics and 
composition of interstellar gas and more detailed studies in¬ 
cluding realistic heating and cooling are required to make 
further progress. 

We hope that this work provides a more general under¬ 
standing of the density variance - Mach number relation in 
the ISM. This is especially true for the warm, atomic part 
of the ISM, where the gas is clearly not isothermal and may 


be approximated with an adiabatic equation of state with 
7 > 1 , which was not covered by previous density variance - 
Mach number relations in the literature. Our new relations 
in this paper attempt to cover this regime and do seem to 
approximately do so, as tested with the set of simulations 
presented here. We hope that the new relations will provide 
useful generalisations of the previous (purely isothermal) re¬ 
lations, which are key ingredients for theoretical models of 
star formation (see e.g., the review by Padoan et al.|[2014|. 
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